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Abstract 


This  paper  investigates  the  asymptotic  distribution  of  the  recently-proposed  contraction  map¬ 
ping  (CM)  method  for  frequency  estimation.  Given  a  finite  sample  composed  of  a  sinusoidal 
signal  in  additive  noise,  the  CM  method  applies  to  the  data  a  parametric  filter  that  matches  its 
parameter  with  the  first-order  autocorrelation  of  the  filtered  noise.  The  CM  estimator  is  defined 
as  the  fixed-point  of  the  parametrized  first-order  sample  autocorrelation  of  the  fdtered  data.  In 
this  paper,  it  is  proved  that  under  appropriate  conditions,  the  CM  estimator  is  asymptotically 
normal  with  a  variance  inversely  related  to  the  signal-to-noise  ratio.  A  useful  example  of  the 
AR(2)  filter  is  discussed  in  detail  to  illustrate  the  performance  of  the  CM  method. 


Abbreviated  Title:  “CM  Method  for  Frequency  Estimation” 

Key  words  and  phrases:  Frequency  estimation,  iterative  filtering,  normality,  parametric 
filter,  spectrum  analysis. 


1  Introduction 


The  contraction  mapping  (CM)  method  is  an  iterative  filtering  approach  for  frequency  estima¬ 
tion  on  the  basis  of  a  noisy  sample  composed  of  a  sinusoidal  signal  in  additive  noise  (see  He  and 
Kedem,  1990a;  Yakowitz,  1991a;  Kedem,  1990;  Li  and  Kedem,  1991).  On  applying  to  the  data 
a  parametric  filter  that  satisfies  the  so-called  fundamental  property,  i.e.,  the  parameter  of  the 
filter  matches  the  first-order  autocorrelation  of  the  filtered  noise,  the  CM  estimator  is  obtained 
by  finding  a  fixed-point  of  the  parametrized  first-order  sample  autocorrelation  of  the  filtered 
data.  In  a  recent  work  by  Li  and  Kedem  (1991),  existence  and  strong  consistency  of  the  CM 
estimator  have  been  proved  based  on  a  general  result  concerning  uniform  strong  consistency  of 
the  parametrized  sample  auto  covariances  of  the  filtered  data.  It  was  shown  (Li  and  Kedem, 
1991)  that  under  some  regularity  conditions  on  the  filter  and  for  sufficiently  large  sample  size, 
the  first-order  sample  autocorrelation  forms  a  contraction  mapping  in  the  vicinity  of  the  true 
frequency,  and  thus  a  unique  fixed-point  can  be  found  with  probability  one  by  using  standard 
iterative  procedures  such  as  the  fixed-point  iteration  (FPI).  Furthermore,  as  the  sample  size 
tends  to  infinity,  the  fixed-point  constitutes  a  strongly  consistent  estimator  of  the  frequency.  In 
this  paper  we  shall  show  that  the  CM  estimator  is  also  asymptotically  normal.  We  shall  also 
provide  a  further  analysis  of  the  AR(2)  filter  considered  by  Li  and  Kedem  (1991),  focusing  on 
its  asymptotic  variance  and  efficiency.  As  a  special  case  of  the  CM  method  using  an  AR(2) 
filter,  Quinn  and  Fernandes  (1988)  showed  that  an  efficiency  as  high  as  that  of  the  nonlinear 
least  squares  (NLS)  estimator  can  be  achieved  in  the  limiting  case  as  the  bandwidth  of  the 
general  AR(2)  filter  shrinks  toward  zero.  This,  together  with  the  results  in  this  paper,  shows 
the  flexibility  and  efficiency  of  the  CM  method  for  frequency  estimation.  A  numerical  example 
is  presented  to  support  the  theoretical  analysis. 

An  entirely  different  proof  of  the  strong  consistency  of  the  CM  method  has  been  given  in 
Kedem  and  Yakowitz  (1991)  using  filtered  zero-crossing  counts  (HOC). 

This  paper  is  arranged  as  follows.  Section  2  briefly  summarizes  the  CM  method  and  some 
asymptotic  properties  obtained  by  Li  and  Kedem  (1991).  In  Section  3,  general  results  are  given 
concerning  the  asymptotic  normality  of  sample  autocovariances  in  the  case  of  a  sinusoidal  signal 
in  additive  colored  noise.  These  results  generalize  those  of  Mackisack  and  Poskitt  (1989)  (see 
also  Stoica  et  ai,  1989)  where  the  asymptotic  normality  was  established  only  for  white  noise. 


1 


The  asymptotic  normality  of  the  CM  estimator  is  presented  in  Section  4  based  on  the  general 
results  obtained  in  Section  3.  Finally,  in  Section  4,  the  AR(2)  filter  is  considered  in  detail  as 
an  important  example  of  the  CM  method.  In  this  section,  the  asymptotic  variance  of  the  CM 
estimator  is  derived  explicitly  for  the  AR(2)  filter,  and  some  numerical  simulations  are  provided 
to  illustrate  its  performance. 

2  The  Contraction  Mapping  Method 

Consider  the  time  series  { yt }  defined  by 

Vt  :=  (3cos(L00t  +  4>)  +  et  (2.1) 

where  /3  >  0  and  cj0  €  (0,  jt)  are  constants,  and  {ct}  is  a  linear  process  of  the  form 

et  =  YJUt-}  (2-2) 

i 

with  (Vbi  <  °°?  and  {6}  being  independent  and  identically  distributed  with  mean  zero  and 
variance  cr|.  The  phase  (f>  of  the  sinusoid  can  be  random  or  constant.  For  convenience,  we  assume 
that  <f>  is  uniformly  distributed  on  (— tt,  7t]  and  independent  of  {£«}.  Under  this  assumption, 
{yt}  becomes  wide-sense  stationary  with  mean  zero.  The  frequency  estimation  problem  is  to 
estimate  the  unknown  frequency  uq  based  on  a  finite  sample  of  {yt}- 

Before  summarizing  the  CM  method,  let  us  introduce  some  useful  notations.  Let  Ca  be 
a  linear  time-invariant  causal  filter  with  real- valued  impulse  response  sequence  {hj(a),j  — 
0, 1, . . .},  where  a  £  [a,  a]  with  a  and  a  being  constants  such  that 

—  1  <  a  <  cos  w0  <  a  <  1. 

Denote  by  H(u;  a)  the  transfer  function  of  Ca,  i.e., 

OO 

H(u-,a):=J2hi(.a)e~iiu- 

j=o 

Let  {yt(a}}  and  (c4(a)},  defined  by 

OO  OO 

Vt(oc)  :=  hj(a)yt-j  and  et(a)  :=  ^  hj(a)et-j, 
j— o  j= o 
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be  the  outputs  of  the  filter  Ca  when  applied  to  {yt}  and  {e*},  respectively.  Define 

.6/^.  E{et+k(a)et{a)} 

^  =  £«(«)} 

to  be  the  fcth-order  autocorrelation  of  {efia)}. 

Suppose  in  the  sequel  that  Ca  is  chosen  such  that  for  any  a  6  [a,  a],  the  following  funda¬ 
mental  property  is  satisfied: 

a  =  p\(a),  (2.3) 

that  is,  the  parameter  a  of  the  filter  matches  the  first-order  autocorrelation  of  the  filtered  noise 
{e4(a)}.  Under  this  assumption,  the  first-order  autocorrelation  of  {yt(a)},  denoted  by  p(a), 
can  be  written  as 


p(a)  —  a*  +  C(a)(a  —  a*) 

where 


(2.4) 


a*  :=  cosuio 


is  the  parameter  to  be  estimated,  and  where 


C(a) := 


1 

1  +  7(a) 


is  the  so-called  contraction  factor  and  7(a)  the  signal-to-noise  ratio  (SNR)  of  {|/f(a)},  i.e., 


7(a)  := 


P2\H(wq',  a)[2 
2  E{e*(a)} 


For  convenience,  let  us  define  G(a )  :=  1  —  C(a)  (known  as  the  gain  factor),  i.e.,  G(ct )  = 
7(a)/{l  +  7(a)},  then  (2.4)  can  be  rewritten  as 


a  —  p(a)  =  (?(a)(a  —  a*).  (2-5) 

It  is  clear  that  if  the  filter  Ca  captures  the  signal  (i.e.,  if  G(a)  >  0,  or,  equivalently,  C(a)  <  1) 
for  all  a  in  a  neighborhood  of  a*,  then  a*  would  be  the  unique  fixed-point  of  the  mapping 
p(a )  in  that  neighborhood.  This  gave  rise  to  the  idea  of  the  CM  method  that  estimates  a*  by 
locating  the  fixed-point  of  an  estimator  of  p(a). 
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The  CM  method  can  be  summarized  as  follows.  For  a  given  sample  {yo,  jq, . . yn-\)  of 
size  n,  let 


t 

yt(tx)  :=  J2hjia)yt-j  (f  =  0,1,..., n-  1). 
j=o 

be  the  filtered  data.  Take  /5(a),  an  estimator  of  p(a),  to  be,  for  instance, 


p(a)  := 


n(q) 

r0(a) 


(2.6) 


where 

n— 1  n— 2 

ro(a)  :=  and  fi(a)  :=  n-1  yt+l(a)yt(a).  (2.7) 

<=o  <=o 

Then,  use  the  so-called  fixed-point  iteration  (FPI) 


4m)  :=  P(4m-1))  (m  =  1, 2,  •  •  •) 


(2.8) 


to  find  a  fixed-point  of  /3(a)  in  the  vicinity  of  a*. 

Assume  now  that  { hj(a )}  satisfies,  in  addition  to  (2.3),  the  following  regularity  conditions 
in  a  closed  subset  A  containing  a*  as  an  interior  point: 

(HI)  There  exist  constants  aj  >  0  such  that 

OO 

y^2jaJ  <  oo  and  \hj(a)\  <  a,j 
j=o 

for  all  j  =  0,1,...  and  all  a  G  A. 

(H2)  h'j(a)  exists  and  there  are  constants  bj  >  0  such  that 

OO 

<  oo  and  |/i'(a)|  <  bj 

3=0 

for  all  j  =  0, 1, . . .  and  all  a  £  A. 

Under  these  conditions,  if  C(a*)  <  1  (i.e,  if  the  signal  is  captured  by  Ca  with  a  =  a*)  and 
the  sample  size  n  is  sufficiently  large,  it  has  been  shown  by  Li  and  Kedem  (1991)  that  almost 
surely  the  mapping  /5(a)  has  a  unique  fixed-point,  denoted  by  d„,  in  a  neighborhood  of  a*,  and 
the  sequence  { 1 )  generated  by  (2.8)  converges  monotonically  to  an  as  m  — *  oo,  provided 
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the  initial  estimate  is  not  too  far  from  an.  Moreover,  the  estimator  an  converges  almost 
surely  to  a*  and  thus  Cjn  defined  by 

ibn  :=  arccosd„  =  arccosp(a„)  (2-9) 

is  strongly  consistent  for  estimating  cjq.  In  the  following,  we  refer  to  both  a„  and  un  as  the 
CM  estimator.  The  main  purpose  of  this  paper  is  to  establish  the  asymptotic  normality  for  an 
and  u»n. 

3  Asymptotic  Normality  of  Sample  Autocovariances 

To  investigcite  the  asymptotic  distribution  of  the  CM  estimator,  we  clearly  need  the  asymptotic 
distribution  of  the  sample  autocovariances  defined  by  (2.7).  For  this  purpose,  let  us  first  consider 
the  following  general  case. 

Let  {hj,j  =  0, 1, . . .}  be  the  real- valued  impulse  response  sequence  of  a  linear  causal  filter 
with  IN  I  <  00 ,  and  let 

OO 

77M  := 

3=0 

be  its  transfer  function.  For  a  given  data  set  {y0,yi,..  .,yn- 1}  obtained  from  (2.1,  define  the 
filtered  data  by 

t 

zt:=Y,hjyt-3  (t  =  0,1,...,  n-  1). 

j=o 

We  would  like  to  establish  the  asymptotic  normality  for  the  sample  autocovariances 

n—j—l 

zt+jzt  (j  =  0,1,..., g)  (3.1) 

<=o 

where  q  is  a  fixed  integer  with  0  <  q  <  n.  For  convenience,  we  assume  that  <f>  is  a  constant.  As 
will  be  seen  later,  this  assumption  can  be  eliminated  without  changing  our  results. 

Note  that  fj  is  defined  for  a  finite  sample.  For  simplicity,  let  us  first  consider  its  counterpart 
fj  defined  on  the  basis  of  an  infinite  sample  by 

71—  1 

■=  n~l  J2  Zl+Cn 

t= o 
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where  zt,  defined  by 


z*  ■=  X  hjyt-h 

3=0 

requires  the  entire  history  of  the  data  {yt}  up  to  time  t.  Using  (2.1),  zt  can  be  rewritten  as 


zt  =  ph  cos(co0t  +<f>h)  +  (t 


where 


Ph  :=  P\H(u0)\,  4>h  :=  <f>+  a,rg{H (lo0)} , 

oo  oo 

Ct  :=  X  =  and  pj  := 

3=0  j  k= 0 


The  following  lemmas  provide  some  asymptotic  properties  of  {fj}. 

Lemma  3.1  Suppose  that  E(£f)  —  ko |  <  oo.  Then 

lim  n  cov( fj,rk)  =  lim  n  T{(fj  -  rj)(rk  -  rk)}  =  Vjk  <  oo 


where 


rj  :=  y  cos(woi )  +  r j ,  rj  :=  E( (t+J (t ) , 


vjk  :=  cos(cuoi)  cos(u;oh)  X  rx  cok(w0t) 

T—  —  OO 
OO 

+  («  -  3)^  +  X  (rTrT+,-fc  +  rT+jrT-k) 


PROOF.  Using  (3.2)  and  the  trigonometric  identity 


COS  U3  cos  A  = 


cos(cu  —  A)  +  cos(tu  +  A) 


rj  can  be  written  as 


X  xJt  +  ^(2n)  1  X  cos(wo(2t  +  j )  +  2<f>h} 
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where 


xjt  ■■=  PhCt+j  cos (u0t  +  fa)  +  PhCt  cos{u>o(i  +  j)  +  fa}  +  Ct+jCt  -  rj.  (3.7) 

Note  that  E(xjt )  =  0  and  that  for  any  j  and  (j>h, 

n— 1 

53  cos{w0(2 1  +  j )  +  2 fa}  =  0(1). 

f=0 

Therefore,  (3.6)  induces 

71—1 

E(fj)  =  rj  +  0(n~1)  and  fj  -  rj  =  n~l  55  X.H  +  Ojn"1).  (3.8) 

t=o 

It  is  readily  seen  from  (3.8)  that 


lim^n  cov(fj,fk)  =  ton  n  E{(fj  -  rj)(fk  -  rk)} 

(n—l  n— 1  \ 

51  xjt,  53 ) 

1=0  1= 0  J 

Moreover,  simple  algebra  shows  that 

(n— 1  n-1  \ 

E  xjtl  5  ]  xkt  I 
1=0  <=0  / 


=  h  +  h  +  Iz 


where 


n—l n—l 


h  :=  Pin"1  Y  55  \rt-s+j-k  cos(w0^  +  fa )  c.os{u0s  +  fa) 

1= 0  s=0 

+  Tt-s+j  cos(u;oi  +  fa)  cos{w0(5  +  k)  +  fa} 

+  As-k  COs{w0(/  +  j)  +  fa}  cos (w05  +  fa) 

+  rj_s  cos{wo(t  +  j)  +  fa}  cos{w0(s  +  k)  +  fa} 

:=  7)  +  T2  +  T3  +  T4, 

n  1  n—l 

h  :=  55W*  “  5  +  h  k)  cos(u0t  +  fa)  +  c(t  -  s,  k)  cos{u>0(t  +  j)  +  fa} 

1=0  5=0 

+  c(s-t  +  k,j)cos(u0s  +  fa)  +  c(s  -  t,j)cos{u)0(s  +  k)  +  fa}], 

in— 1  n—l  \ 

1  cov  (  55  C*+y  Ct ,  53  Ct+kCt ) » 


I: 


n 


\t=0 


t- 0 


/ 


and 


c(jik)  7?(Ct+jCt+fcCi) 
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is  the  third-order  cumulant  of  {Ct}*  Using  (3.5)  and  the  substitution  r  :=  t  —  s,  '1\  can  be 
written  as 


T,=&  V  rC 

1  cy  VT+j-k 

|r|<n— 1 


1  -  — )  cos(cj0t)  +  -  53  cos{o;o(2t  -  r)  +  2  fa} 

nJ  ntn 


where  D  :=  {t  :  max(0,r)  <  t  <  min(n  —  l,r  +  n  —  1)}.  Since  lr$l  <  °°  and 


n  1  53  cos{cjo(2t  —  r)  +  2fih}  — »  0  as  n  — *■  oo 

teD 


for  any  r,  bounded  convergence  theorem  gives 


Pi 

2 


2  oo  o' 

C  _  P i 


2  oo 


T1  ^  X]  rr+j-/t  c°s(^or)  =  X]  ^COs{u>0(t- j  +  fc)}. 


Similarly,  we  obtain 


T3 

t4 


Pi 


2  oo 


r$  cos{o;o(r  -  j  -  fc)} 


Pi 


T  —  —  00 

2  o° 


53  cos{w0(r  +  j  +  k)} 


T—  —  00 

2  oo 


%  Yj  rrCos{a >o(T  +  j-k)}. 


Combining  these  together  yields 

OO 

h  ->■  2Plcos(u>0j)cos(L00k)  53  tt  cos(wo'r)- 


In  addition,  since  |c(r  +  j,k) |  <  oo  for  any  fixed  j  and  k,  an  analogous  argument  leads  to 
I2  — >  0.  Finally,  Proposition  7.3.1  (Brockwell  and  Davis,  1987)  gives 

OO 

J3 -►  (/c  -  3)r£r£  +  53  (rir(T+j-k  +  4+Yr-k)  ■ 

r— — 00 

The  assertion  is  thus  proved  upon  combining  these  results.  <)> 


Lemma  3.2  Under  the  same  assumption  as  Lemma  3.1,  {nl!2[fj—rj),j  =  0, . .  .,q}  are  asymp¬ 
totically  jointly  normal  with  mean  zero  and  covariance  matrix  V  ( Vjk)9j  k=o >  where  vjk  is 
defined  by  (3.4). 
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Proof.  Using  the  Cramer- Wold  device  (Brockwell  and  Davis,  1987,  Proposition  6.3.1),  we  only 
need  to  show  that  for  any  {Ay,  j  =  0, . .  {0}, 

9 

n1/2  ^  A j(rj  -  rj )  — »■  N(0,  v)  in  distribution, 
j= o 

where  i >  :=J2Y1  AjA kVjk  >  0.  From  (3.8),  it  suffices  to  show  that 

71—1 

n-1'2  ^  xt  — »  N(0,  v)  in  distribution,  (3-9) 

t= o 

where  xt  :==  J2  ^jxjt->  and  Xjt  is  defined  by  (3.7).  If  the  sequence  {/r,}  given  by  (3.3)  has  a 
finite  length,  i.e.,  if  /Uy  =  0  for  all  |j|  >  m  with  m  being  a  positive  integer,  then,  after  some 
straightforward  manipulations  similar  to  the  proof  of  Theorem  6.7  (Hall  and  Heyde,  1980,  pp. 
192),  we  can  write 

n-i/2  Y^xt  =  n~1/2Y2Mt  +  oP(l) 

t= o  t= 1 

where 

2  m+q  q 

Mt  :=  Bt£t  +  ^  AT(£t£t-T  —  At  :=  A ?A?T, 

r=o  i=o 

9 

-B'f  :=  Y'^/bfc/^[c°s{wo(t  +  fc  +  i)  +  +  cos(wo(t  +  A;  -  j)  +  <^J], 

j= 0  k 

Aj0  :=  rj/cr2,  and  Ajr  :=  (rJ+T  +  rj_T)/a|,  for  r  ^  0. 

It  is  easy  to  see  that  {M<}  are  martingale  differences  with  respect  to  the  a-fields  J~t  generated 
by  {£j,i  <  t}  f°r  t  >  0.  Direct  computations,  followed  by  an  application  of  the  law  of  large 
numbers,  yield 

n 

n~xV2  :=  n_1  E(M2\Jrt-i)  —*  v  in  probability. 

t= l 

It  can  also  be  shown  by  direct  computations  that 

n_14  :=  n~xE{V2)  -*•  »■  (3.10) 

Therefore,  l/r2  /.s2  — *  1  in  probability.  Furthermore,  (3.10)  imphes  that  the  Lindeberg  condition 
(Brown,  1971,  eq.  2)  is  equivalent  to 

n 

n-'YsHMflilMtl  >e*)}->0, 

t=i 
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where  /(•)  is  the  indicator  function.  Hence  it  suffices  to  verify  that 


E{M?I(\Mt\  >  £*)}  -  0 


for  any  e  >  0.  To  do  so,  we  note  that  \Bt\  <  B  for  some  B  >  0  and  all  t.  Therefore, 

\Mt\  <  Ut  :=  5|6|  + 

Moreover,  (3.10)  implies  that  st  >  et1/2  for  small  e  and  large  t.  Combining  these  results  gives 
E{M2I(\Mt\  >  £5<)}  <  E{U2I{Ut  >  eh1!2)}  =  E{U2I(U0  >  eh1!2)}  -  0, 


2  m-\-q 

T.  AT((it£t-T  — 

t= 0 


where  the  second  expression  is  due  to  the  stationarity  of  {Ut}.  Applying  Theorem  2  (Brown, 
1971)  proves  (3.9). 

If  {/,„}  is  of  infinite  length,  (3.9)  can  also  be  verified  by  following  the  proof  of  Proposition 
7.3.3  (Brockwell  and  Davis,  1987).  In  fact,  for  any  m  >  0,  let  us  define 


/■m 

St 


~m 

r3 


y  HjZt-j,  Z?  :=  fh  cos (w0f  +  4>h)  +  cr> 

-m 


j——m 

n—1 

—1 
n 

t= 0 


If  so  doing,  we  obtain,  for  any  fixed  m, 

j= 0 


and 


n—1 

S™  :=  n"1/2  y  x™  -»  Sm  ~  N(0,um) 

t= o 

in  distribution  as  u  — >  oo,  where  x™  :=  Y1  ^jx]t  and  x]t  is  given  by  (3.7)  with  in  place  of  (*, 
and  where  vm  :=  X)  Y.  and  vj‘k  is  defined  by  (3.4)  with  autocovariances  of  {£*}  replaced 

by  autocovariances  of  {C™}.  Since  vm  — *■  v  as  m  —»•  oo,  then  Sm  — ►  N(0,u)  in  distribution. 
Moreover,  an  analogous  argument  as  in  the  proof  of  Proposition  7.3.3  leads  to 


lim  limsup  x>v{nll2\r?  —  rf  —  f\.  +  r,-|  >  e)  =  0 

m-*  OO  J  1 

for  any  e  >  0.  The  proof  is  then  completed  by  applying  Proposition  6.3.9  (Brockwell  and  Davis, 
1987).  <> 
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Remark  3.1  The  asymptotic  normality  has  been  proved  by  Mackisack  and  Poskitt  (1989)  for 
the  simplest  case  where  {£<}  are  independent  and  identically  distributed  (see  also  Stoica  el  al. , 
1989).  Therefore,  Lemma  3.2  is  an  extension  of  their  results  to  the  case  of  colored  noise.  § 

With  the  help  of  the  above  lemmas,  we  are  now  able  to  state  the  central  limit  theorem  for 
the  sample  autocovariances  {fj}  given  by  (3.1). 

Theorem  3.1  Suppose  that  E(£f)  =  ko |  <  oo  and  that 

OO 

H iN  <  oo- 

i=o 

Then,  —  rj'),  j  =  0, . . . ,  q}  are  asymptotically  jointly  normal  with  mean  zero  and  covari¬ 

ance  matrix  V. 


Proof.  By  Lemma  3.2,  it  suffices  to  show  that  n1/,2(rT  —  rT )  — >■  0  in  probability  as  n  oo.  To 
do  so,  we  first  claim  that 

nl!2(fT  —  fT)  — *•  0  in  probability,  (3.11) 


where 


n— 1 


rV  :=  n  1  Y,  Zt+rZt- 


t= 0 


In  fact,  simple  algebra  gives  zt  =  zt  +  zt  and 


n—  1 

n1/2(fT-rT )  =  n~1/2  Y^(zt+Tzt  +  zt+Tzt  +  Zt+rZt) 

t~  o 

•=  h  +  I2  +  I3 


where 

OO 

Zt  ^2  hjyt-j. 
j=t+i 

Using  (2.1),  Ii  can  be  written  as 

71—1  OO  t~\~T 

h  -  n~1/2J2  Y2  Y2hjhk\fl2  c°s{LV0(t  -  j)  +  </>}  cos{u0(t  +  r  -  k)  +  <j>} 

<= 0  j=t+ 1  jt=0 

+  /3et-j  cos{Lo0(t  +  t  -  k)  +  <fi} 

+  fdet+r-k  cos{w0(f  -  j)  +  <t>} 

+  Q— jO+T— fc] 

:=  T1  +  T2  +  T3  +  T4. 
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Clearly,  we  have 


71—1  OO  OO 

| Ti|  <  n  Y  Y  N  <  n~1/2/?2i7  £ j'N  ^  0  as  n  -  oo, 

t=o  j=t+ 1  j= o 

where  H  '■=  J2\hj\-  Similarly, 

CO 

T(|T2|)  <  rT1/2  Edfoiy^/Eil^l  -  0  as  n  -  oo, 

i=o 

and  thus  T2  — >  0  in  probability  by  applying  Chebyshev’s  inequality.  In  the  same  way,  we  can 
show  that  T3  and  T4  vanish  in  probability  as  well.  Combining  these  results  gives  I\  — *  0  in 
probability.  The  same  conclusion  can  also  be  made  about  /2  and  /3  by  an  analogous  argument. 
(3.11)  is  thus  proved.  Finally,  we  claim  that  n1/2(rT  —  fT)  — ►  0  in  probability.  This  can  be 
easily  established  upon  noting  that 

n— 1  t  t-\-r 

n1/2  E(\fT  -  fT |)  <  n”1/2  E  E  E  I^AI  T(|2/<+T_fet/t_j|)  <  nT^rKH2 

t—n-T  j= 0  A;=0 

where  the  last  inequality  is  due  to  the  fact  that 

E(\ytys\)  <  K  :=  /32  +  2/?f?(|€o|)  +  £(*o)- 

The  theorem  is  thus  proved.  () 

Remark  3.2  Theorem  3.1  remains  valid  if  4>  is  random  and  independent  of  {£t}.  In  fact, 
by  conditioning  on  <f>,  Theorem  3.1  can  be  proved  in  the  same  way.  Since  the  asymptotic 
conditional  distribution  given  (f>  does  not  depend  on  <f>,  the  same  result  holds  for  the  asymptotic 
unconditional  distribution  as  well.  Q 

Remark  3.3  This  theorem  can  be  easily  generalized  to  the  case  of  multiple  sinusoids  in  additive 
noise.  In  fact,  suppose  that  {yf}  is  given  by 

9-1 

yt  ■=  E  P"  cos(w^  +  4u)  + 

i/=0 

where  q  is  a  postive  interger,  /?„  and  are  constants  with  /?„  >  0  and  0  <  <  •  •  •  <  uq  <  7r, 

and  (f>„  6  (— ir,  7r]  can  be  either  constants  or  random  and  independent  of  the  noise  {t,*}  given  by 

(2.2).  Then,  Theorem  3.1  still  holds  with  Vjk  defined  by 

q— 1  (  00 

vjk  :=  2  Y  i  Plh  cos(uv))  cos^fc)  Y  ri  cos(w1/r  ) 

i/—0  K  T——00 

00 

+  (fC-3)r|f<+  Y  {^r+j-k  +  rr+/r-k) 

T-  —  CO 
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where  f3uh  :=  /3\H(u„)\.  0 

By  a  somewhat  different  approach,  Yakowitz  (1991b)  have  given  a  result  similar  to  Theo¬ 
rem  3.1.  But  in  his  work,  only  a  lower  bound  of  the  asymptotic  variance  is  provided  and  the 
result  is  restricted  to  FIR  filters. 

As  a  result  of  Theorem  3.1,  we  can  also  establish  the  asymptotic  normality  for  the  sample 
autocorrelation 

Pi  '■-=■  Tj/ro 

as  follows. 

Corollary  3.1  Suppose  that  conditions  in  Theorem  3.1  are  satisfied.  Then,  nl/2([>j  —  pj )  is 
asymptotically  normal  with  mean  zero  and  variance  Vj,  where  pj  :=  rjjr o  and 

vi  ■=  (p]v oo  ~  2 pjVQj  +  Vjj)/rl.  (3.12) 


Proof.  The  assertion  follows  from  Theorem  3.1  by  applying  Proposition  6.4.3  (Brockwell  and 
Davis,  1987).  <) 

4  Asymptotic  Normality  of  the  CM  Estimator 

Applying  the  general  results  obtained  in  the  proceeding  section,  we  would  like  in  this  section 
to  establish  the  asymptotic  normality  for  the  CM  estimators  an  and  uin  defined  in  Section  2. 
For  this  purpose,  we  first  have 

Lemma  4.1  Let  {yt}  be  defined  by  (2.1)-(2.2)  and  assume  that  <  oo.  Suppose  in  addition 

that 

OO 

^2j\ hj(a*)\  <  oo. 

3=0 

Then,  n1/2{p(a*)  --  p(a*)}  is  asymptotically  normal  with  mean  zero  and  variance  o 2p,  where 
p{ct*)  is  given  by  (2.6)  with  a  —  a* ,  and  o2p  :=  vpC2(a*)  with 

OO 

”,■=  E  {(l  +  2o*2)(4(«*))2-4a*4(a-)4+I(a-)  +  4+1(aK_,(«)}.  (4.1) 
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Proof.  Applying  Corollary  3.1  with  hj  hj(a*)  and  (t  ■=  e*(a*),  we  are  able  to  obtain 


n1/2{p(a*)  -  p(a*)}  ->  N(0,  vj)  in  distribution, 

where  rq  is  defined  by  (3.12)  with  j  =  1.  Note  that  in  (3.12),  r0  =  a 2/C(a*)  where  a\  := 
E((f)  =  rg.  Note  also  that  p(a*)  =  p\(a*)  =  a*  =  cosu>0.  Therefore,  (3.4)  gives 

OO 

W =  S  +  2  J2  ( Pk(a *))2 

k— — oo 
oo 

W®c  =  u*S  +  2  Pk(a*)Pk+ i(tt*) 

k=—oo 

oo 

Vii/v*  =  a*2S+  {(/>fc(“*))2  +  Pfc+i(«*K-i(«*)} 

k— — oo 

where  S  :=  (2/32/<72)  p|(a*)  cos(u?oA;)  -f  k  —  3.  Simple  algebra  yields  v\  =  <r2  <$> 

Remark  4.1  Because  of  asymptotic  equivalence,  there  are  many  other  ways  of  estimating  p(a) 
without  changing  the  results  in  Lemma  4.1.  In  fact,  p(a )  can  be  taken,  for  example,  as  the 
least  squares  estimator  that  minimizes 

£{&(«)  ~  PVt- 1(«)}2  +  YliVt- 2(a)  -  pyt-i(a)}2. 

t= 1  t= 1 

In  this  case, 

n— 1 

J2  yt-i(a){yt{a)  +  yt-2{ot)} 

p(a)  =  ^ - — - .  (4.2) 

2Eti(«) 

t=i 

Some  other  estimators  of  p(a)  can  be  found  in  Li  and  Kedem  (1991).  <0> 

Using  this  lemma,  we  are  able  to  obtain  the  following  central  limit  theorem  for  the  CM 
estimator: 


Theorem  4.1  Let  {yf\  be  defined  by  (2.1)-(2.2)  and  suppose  that  E(£f)  <  00.  Assume  further 
that  (HI),  (H2)  are  satisfied  and  that  G(a )  >  g  >  0  for  all  a  in  a  neighborhood  of  a*.  Then, 
nll2(an  —  o:*)  and  nl!2(un  —  u;0)  are  asymptotically  normal  with  mean  zero  and  variances 


72(a*) 


and 


= 


72(cosu;o)  sin2o;o  ’ 


(4.3) 


respectively. 
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Proof.  Let  en(a)  :=  /5(a)  —  p(a).  Since  an  is  a  fixed-point  of  p(a),  then  p(an)  =  an. 
Combining  this  with  (2.5)  yields 

®n(®n.)  —  G(an)(an  a  ). 

Under  assumption  (H2),  en(a )  is  differentiable  with  probability  one.  Therefore,  we  have  the 
Taylor  expansion  G(an)(dn  -  a*)  =  en(a*)  +  e'n(an)(an  -  a*),  or, 

{G(d„)  -  e'n(an)}(an  -  a*)  =  e„(a*),  (4.4) 

where  an  lies  between  a*  and  an.  According  to  Theorem  5.1  and  Corollary  5.1  (Li  and  Kedem, 
1991),  we  have,  almost  surely,  an  — >  a*,  and  e'n(a)  =  p'(a)  —  p'(a)  —>  0  uniformly  in  a 
neighborhood  of  a*.  Therefore,  by  the  continuity  of  G(a),  we  obtain 

G(an )  —  e'n(an )  — >•  G(a*)  in  probability. 

Moreover,  Lemma  4.1  guarantees  that  n1/2e„(a*)  is  asymptotically  normal  with  mean  zero  and 

variance  a2p.  The  asymptotic  normality  of  an  with  a\  :=  a2p/G2(a*)  is  thus  proved  by  combining 

all  these  results  and  applying  Slutsky’s  theorem  to  (4.4).  From  (4.1)  and  the  definitions  of  G(a ) 

and  C(a),  we  obtain  a\  =  vpj^2{a*).  Finally,  since  cb„  =  arccosdn  and 

d  arccos  a  —  1 

da  a=a‘  |sinwo|’ 

the  asymptotic  normality  of  un  follows  from  Proposition  6.4.1  (Brockwell  and  Davis,  1987).  0 
Remark  4.2  The  formula  (4.3)  shows  that  in  order  to  reduce  the  asymptotic  variances,  the 
filter  must  be  chosen  so  as  to  enhance  the  sinusoidal  signal  when  a  takes  on  the  true  value  a* . 

❖ 

5  The  AR(2)  Filter 

Consider  the  following  AR(2)  filter  defined  recursively  by 

yt(a)  +  0(a)Wt-i(°O  +  V2yt-2(a)  =  yt  (t  =  0, 1, . . .,  n  -  1),  (5.1) 

where  0  <  p  <  1,  and 

1  -I-  Tt2 

9(a)  := - -G-a.  (5.2) 
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For  any  |a|  <  277/ (1  +  if),  since  |0(a)|  <  2,  we  can  alway  write 

0(a)  —  — 2cos{u>(a)}  (5-3) 

for  some  u(a)  €  [0,7r].  It  turns  out  that  ?7exp{±iu>(a)}  are  poles  of  the  AR(2)  filter  defined  by 
(5.1).  Moreover,  it  can  be  verified  (Li  and  Kedem,  1991)  that  the  fundamental  property  (2.3) 
is  satisfied  by  this  filter  whenever  {et}  is  white  noise,  i.e.,  whenever  e*  are  independent  and 
identically  distributed  with  mean  zero  and  variance  of .  Therefore,  in  this  case,  the  filter  (5.1) 
can  be  applied  to  estimate  the  frequency  oj0  by  using  the  FPI  procedure  (2.8)  and  by  taking 

:=  arccos/5(d^m-1))  =  arccosd^”1^  ( m  =  1,2, . . .).  (5-4) 


Here  p(a)  can  be  defined  by  either  (2.6)  or  (4.2),  and  yt(a )  is  given  by  (5.1)  with  zero  initial 
values.  This  procedure  is  clearly  an  iterative  filtering  algorithm.  In  fact,  at  the  mth  iteration, 
we  first  filter  the  sample  { yt }  by  (5.1)  with  a  —  and  then  calculate  /5(dlm  ^)  on  the 

basis  of  the  output  using  (2.6)  or  (4.2).  Finally,  we  obtain  by  (5.4)  as  an  estimate  of  u0 
and  a^n'>  by  (2.8)  for  the  next  iteration. 

It  has  also  been  verified  (Li  and  Kedem,  1991)  that  the  filter  (5.1)  satisfies  (III)  and  (H2) 
with  C(a*)  <  1.  Therefore,  it  is  guaranteed  with  probability  one  for  large  n  that  the  mapping 
/5(a)  has  a  unique  fixed-point  an  in  the  vicinity  of  a*,  and  that  defined  by  (5.4)  converges 
monotonicadly  to  un  —  arccosan  as  m  — »  00.  Furthermore,  Cjn  is  strongly  consistent  for 
estimating  uq  and,  according  to  Theorem  4.1,  n1//2(o>„  —  u>o)  is  asymptotically  normal  with 
mean  zero  and  variance  of  defined  by  (4.3). 

To  obtain  explicit  expressions  for  erf  and  of, ,  we  first  note  that  for  the  AR(2)  filter  (5.1), 


*\|2 


\H(u0-,a*)\ 


1 


(1  —  rj2)2  sin2  uq 

Using  the  formulas  given  by  He  and  Kedem  (1990b)  for  AR(2)  filters,  we  obtain 


sin2M«)} 


J2  sin2{o;(a)(y  +  1  )}t]2j, 


and 


1  —  772 

p'^(«)  =  cos{n>(a)A;}?7^  4-  y-jj — —  cot{w(a)}  sin{ta(a)fc}?7^ 


(5.5) 
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for  k  >  0,  where  u(a)  is  defined  by  (5.3).  It  follows  by  straightforward  computations  that 


E{c2t(a)}  = 


a 


2 

£ 


(  1  -  774)(1  -a2)’ 


and,  therefore,  we  can  write 


7(«*)  = 


1  +  rf 

T^1 


where  7  :=  /?2/(2 a2)  is  the  signal-to-noise  ratio  (SNR)  of  {yt}-  It  is  readily  seen  that  7(0*), 
the  signal-to-noise  ratio  of  the  filtered  data  at  a*,  does  not  depend  on  a*.  This  means  that 
the  sinusoidal  signal  is  equally  enhanced  by  the  AR(2)  filter  (5.1)  for  any  u>0  G  (0,7r)  satisfying 
|  cos^ol  <  277/(1  +  T)2). 

In  addition,  direct  but  lengthy  calculation  of  (4.1)  using  (5.5)  yields 


l-T}2 

'<' =  TT?(1  “  0  >• 


Therefore,  from  (4.3),  we  obtain 


1  —  i]2  y  1  —  a 

I  +  772, 


*2 


T 


and 


(5.6) 


Clearly,  in  order  to  achieve  a  small  asymptotic  variance,  77  should  be  chosen  as  close  to  1  as 
possible. 

Given  (5.6),  a  natural  question  is  what  we  can  learn  from  these  formulas  about  the  mean- 
square  error  (MSE)  E(un  —ujo)2  of  the  CM  estimator  u?n.  Since  a  closed-form  expression  of  the 
MSE  is  difficult,  if  not  impossible,  to  obtain,  we  would  like  to  investigate  this  problem  based 
on  a  numerical  experiment.  First  of  all,  it  should  be  noted  from  (5.6)  that  a2  tends  to  zero  as 
77  — >  1.  Therefore,  it  is  reasonable  to  believe  that  in  order  for  the  asymptotic  variance  a2  jn  of 
un  to  be  a  good  approximation  of  the  MSE,  a  larger  sample  size  n  is  required  when  77  becomes 
closer  to  1,  as  vindicated  by  the  simulation  results  shown  in  Table  1. 

For  various  values  of  77  and  n,  Table  1  presents  the  estimated  MSE  of  the  norm alized  CM 
estimator  d>n/7r  on  the  basis  of  100  independent  realizations.  Each  realization  was  generated 
from  the  model  (2.1)-(2.2)  with  {e*}  being  Gaussian  white  noise,  and  with  lo0  =  0.42x,  =  O-lx, 

and  7=1  (i.e.,  SNR  =  0  dB).  Figure  1  gives  the  corresponding  curves  obtained  from  Table  1. 
To  produce  the  CM  estimates,  the  FPI  procedure  (2.8)  was  used  in  connection  with  p(a)  given 
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Table  1:  MSE  of  CM  Estimator  by  FPI  with  AR(2)  Filter  (SNR=C 


0.85 

7.25  XlO-6 
4.24x10-® 
1.71X10-6 
1.41x10-6 
8.97  x  10“7 
8.47  XlO"7 
8.29  X  10-7 
6.05  XlO-7 
5.43  XlO"7 
4.71  XlO”7 
3.58  XlO"  7 
3.85  XlO- 7 

3.49  X 10“  7 

3.26  XlO-7 
2.31  Xl0~7 
2.82  XlO-7 
2.33  Xl0“7 

2.49  XlO"  7 
2.21  XlO-7 
2.23  xlO“7 


3.06x10-® 
1.17xl0-e 
5.35  XlO-7 
3.91  XlO-7 
2.82  XlO-7 
2.34  xlO~7 
2.47x10— 7 
1.67  X 10-7 
1.67  X  10-7 
1.30  XlO"7 
1.05  XlO-7 
1.07  X  10_ 7 
9.66X10-8 
9.02X10-8 
7.54x10-8 
7.81  X  10-8 
6.83X10-8 
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Figure  1:  —  Log  plot  of  MSE  of  the  CM  estimator  for  ( from  the  bottom )  r]=0.85,  0.90,  0.95, 

0.98,  0.99,  and  0.999.  The  dashed  lines  ( - )  are  asymptotic  variances  of  the  CM  estimator 

for  ( from  the  bottom )  p—0.85,  0.90,  0.95,  0.98,  and  0.99.  The  solid  line  is  the  asymptotic 
variance  of  the  NLS. 

by  (4.2).  For  p  =  0.85,  the  initial  value  was  taken  to  be  =  cos(0.67r)  (i.e.,  u>i°^  =  0.67 r).  The 
FPI  was  terminated  according  to  the  stopping  criterion:  m  <  20  or  \u>^  —  wim  ^|  <  7T  X  10-5, 
and  the  resulting  estimate  was  used  to  initiate  the  FPI  for  the  next  (larger)  value  of  p.  The 
corresponding  values  of  <y2/(ir2n)  are  shown  in  every  second  line. 

Two  conclusions  can  be  drawn  immediately  from  these  results,  (a)  Given  a  fixed  p,  the 
closer  it  is  to  1,  the  larger  the  sample  size  n  is  required  for  er^/n  to  be  a  good  approximation 
of  E {Cjn  —  wo)2-  For  instance,  if  p  =  0.85,  good  approximations  are  given  for  n  >  300;  while 
for  p  =  0.95,  we  need  n  >  1700.  (b)  Given  a  fixed  sample  size  n,  the  MSE  can  be  reduced  by 
increasing  p.  In  the  limiting  case  as  p  — >  1,  the  CM  estimator  achieves  the  same  efficiency  as 
does  the  nonlinear  least  squares  (NLS)  estimator.  Indeed,  as  shown  by  Quinn  and  Fernandes 
(1988),  in  the  case  of  p  =  1  we  have 

n3/2(w„  —  wo)  — »  N(0,<7^r£S)  in  distribution, 

where  :=  12/p.  Therefore,  for  fixed  n,  the  asymptotic  variance  offLS/n3  of  the  NLS 

estimator  can  be  utilized  as  an  approximation  of  the  MSE,  if  p  is  very  close  to  1. 
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Finally,  it  should  be  pointed  out  that  the  Quinn-Fernandes  procedure  (a  special  case  of 
the  CM  method  using  the  FPI  and  the  AR(2)  filter  with  /3(a)  given  by  (4.2)  and  i)  =  1) 
requires  the  initial  estimate  to  be  of  accuracy  o(n-1/2)  (see  Quinn  and  Fernandes,  1988),  and 
therefore  they  suggested  using  some  other  methods  to  initiate  the  procedure.  On  the  other 
hand,  the  CM  method  with  r/  <  1  requires  less  stringent  initial  estimates.  (In  Table  1,  Au  := 
—  cj0|  =  0.18?r.)  In  fact,  the  CM  estimates  obtained  with  a  smaller  rj  can  be  utilized  as 
initial  estimates  for  a  larger  rj.  Eventually,  as  rj  — >  1,  one  can  still  achieve  the  same  efficiency 
of  the  NLS  estimator.  In  conclusion,  the  CM  method  provides  a  flexible  and  efficient  approach 
of  frequency  estimation. 
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